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EXPLICIT  DIFFERENCE  SCHEMES  FOR  WAVE 
PROPAGATION  AND  IMPACT  PROBLEMS* 


Joseph  E.  Flaherty 
Department  of  Mathematical  Sciences 
Rensselaer  Polytechnic  Institute 
Troy,  NY  12181 

and 


U.S.  Army  Armament  Research  and  Development  Command 
Large  Caliber  Weapon  Systems  Laboratory 
Benet  Weapons  Laboratory 
Watervliet,  NY  12189 


ABSTRACT.  Explicit  finite  difference  and  finite  element  schemes  are 
constructed  to  solve  wave  propagation,  shock,  and  Impact  problems.  The 
schemes  rely  on  exponential  functions  and  the  solution  of  linearized  Rlemann 
problems  in  order  to  reduce  the  effects  of  numerical  dispersion  and  diffusion. 
The  relationship  of  the  new  schemes  to  existing  explicit  schemes  is  analyzed 
and  numerical  results  and  comparisons  are  presented  for  several  examples. 

I.  INTRODUCTION.  Exponentially  fitted  and/or  weighted  finite  difference 
19, 111,  finite  element  [3,9,10],  and  collocation  [5]  schemes  have  become 
popular  and  effective  methods  of  solving  steady  convection-diffusion  problems. 
They  avoid  the  spurious  mesh  oscillations  that  are  found  near  boundary  and 
shock  layers  when  centered  schemes  are  used  at  high  cell  Reynolds  or  Peclet 
numbers  and  they  reduce  the  effects  of  numerical  diffusion  that  are  associated 
with  classical  upwind  difference  schemes. 

We  seek  to  extend  exponential  methods  to  transient  problems  and  as  a 
first  step  we  consider  one-dimensional  scalar  initial  value  problems  of  the 
form 


ut  +  f(u)x  -  eu**  ,  t  >  0  ,  |x|  <  • 

u(x,0)  -  u°(x)  ,  |x|  <  •  (1) 

where  0  <  e  «  1  is  either  a  real  or  an  artificial  viscosity  parameter  and  the 
x  and  t  subscripts  denote  partial  differentiation. 

Our  primary  motivation  for  studying  exponential  schemes  is  a  desire  to 
develop  improved  numerical  methods  for  elastic-plastic  Impact  problems  in 
solids  and  blast  problems  in  gases. 


*This  research  was  partially  sponsored  by  the  U.S.  Air  Force  Office  of 
Scientific  Research,  Air  Force  Systems  Command,  DSAF,  under  Grant  Number 
AFOSR  80-0192.  The  United  States  Government  is  authorized  to  reproduce  and 
distribute  reprints  for  government  purposes  notwithstanding  any  copyright 
notation  thereon. 
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In  this  paper,  we  confine  our  attention  to  explicit  difference 
approximations  of  (1)  having  the  fora 

1 

Un+lj  -  D°j  -  -  11(1  +  "  ^J-l)  +  (1  “  *nJ+l/2)(fnJ+l  ~  f°j)l 

| 

+  —  (U°j-i  -  20°j  +  Dnj+i)  ,  n  >  0  ,  |J|  <  - 

0°j  -  U°( J  Ax)  ,  |j|  <  -  (2) 

where  Ax  and  At  denote  the  uniform  spatial  and  temporal  grid  apaclngs, 
respectively,  U°j  Is  the  numerical  approximation  of  u(jAx,nAt), 

fnj  f(Onj)  ,  X  -  At/ Ax  (3,4) 

and  znj+i/2  are  upwind  weighting  factors. 

Many  popular  difference  schemes  have  the  form  of  (2)  and  some  of  these 
are  discussed  and  compared  In  Section  II.  We  also  Introduce  an  exponential 
scheme  that  Is  based  on  determining  *nj±i/2  *°  that  U°j  is  the  exact  solution 
of  the  linearized  steady  equation 

cu*  “  cuxx  (5) 

when  c  -  f'(u)  is  a  constant.  We  call  this  method  the  linearized  steady 
exponential  (LSE)  method  and  It  is  the  simplest  extension  of  the  exponentially 
fitted  and  weighted  schemes  (3,9,10,11]  to  transient  problems.  The  scheme 
gives  Improved  accuracy  for  steady  shock  problems,  but  offerB  little  Improve¬ 
ment  over  classical  upwind  differencing  for  moving  shock  problems. 

In  Section  III  we  develop  an  exponential  scheme  that  Is  based  on  the 
,  exact  solution  of  a  linearized  transient  equation  (1)  that  is  subject  to 

piecewise  constant  Initial  data,  l.e.,  a  linearized  Riemann  problem  for  (1). 

We  call  this  method  the  linearized  transient  exponential  (LTE)  method  and, 
like  other  methods  based  on  the  solutions  of  Riemann  problems  [1,2,4,6,8,13, 
t  15,16],  It  sharply  resolves  boundary  and  shock  layers  without  added  diffusion 

or  spurious  oscillations.  As  t  ♦  0  the  LTE  method  becomes  formally  equivalent 
l  ,  to  Roe's  method  [15,16]  for  hyperbolic  systems  of  conservation  laws,  van  Leer 

£  [17]  has  noted  that  Roe's  method  treats  an  expansion  fan  as  a  so  called 

*  "expansion  shock"  (cf.  Figure  6)  and,  unfortunately,  our  LTE  method  also  has 
this  disturbing  property  even  when  c  is  nonzero,  but  small. 

*  . 

In  Section  IV  we  present  some  preliminary  results  for  vector  systems  of 
equations  and  in  Section  V  we  discuss  our  results  and  indicate  some  directions 
for  future  work. 

‘  II.  THE  LINEARIZED  STEADY  EXPONENTIAL  (LSE)  METHOD.  The  LSE  method  is 

obtained  from  (2)  by  selecting  z1^,  k  ”  J ±1/2 ,  as 

•  coth  p\/2  -  2/p\  (6) 

I 

I 
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L  j 


0 


where  P°k  !•  the  cell  Reynold*  number 


Ax  <^+1/2  ♦  cVl/2 

p'ht  -  -  ( - : - )  (7) 

t  2 

and 

f'O^j)  (8) 

As  previously  noted,  the  LSE  method  will  give  a  pointwise  exact  steady  state 

solution  of  (1)  when  cnj  is  a  constant.  This  or  similar  schemes  have  been 
used  by  several  investigator*  [3,9,10,11]  for  steady  singularly-perturbed 
problems  and  herein  we  try  to  apply  it  to  transient  problems. 

We  first  consider  a  linear  stability  analysis  of  the  difference  scheme 
(2)  by  letting  f(u)  ■  cu  where  c  is  a  constant.  We  also  let  p  and  z  denote 
the  constant  values  of  and  z\,  respectively,  a  denote  the  Courant  number 

a  -  cAt/Ax  (9) 

and  0  denote  the  dissipation  parameter 

0  -  a(z  +  2/p)  (10) 

In  this  case,  equation  (2)  can  be  written  as 
..  1  1 

Un+ij  -  unj  -  -  a(Dnj+1  -  U^)  +  -  0(Unj+1  *  2Unj  +  Unj-i)  (11) 

Several  popular  difference  schemes  have  the  form  of  (11)  for  different 
values  of  0  (or  z)  and  some  of  these  are  listed  in  Table  1.  All  of  these 
schemes  are  first  order  accurate  in  time,  except  the  Lax-Wendroff  scheme  which 
is  second  order. 

A  von  Neumann  analysis  (cf.  Richtmver  and  Horton  [14])  shows  that 
equation  (11)  is  stable  in  the  region  a*  <  0,  0  <  0  <  1.  This  region  is  shown 
shaded  for  o  >  0  in  Figure  1.  Curves  corresponding  to  the  methods  in  Table  1 
are  also  shown.  We  see  that  the  LSE  method  slightly  improves  upon  the 
stability  and  accuracy  properties  of  upwind  differencing  and  that  centered 
differencing  and  the  Lax-Friedrichs  scheme  are  outside  of  the  stability  region 
for  most  values  of  a  and  0. 

Example  1:  We  compare  the  methods  in  Table  1  on  the  constant  coefficient 
initial-boundary  value  problem 


ut  +  ux  "  Cuxx  »  t  >  0  ,  0  <  x  <  1 

u(x,0)  -  0  ,  0  <  x  <  1 

u(0,t)  -  1  ,  u( 1 , t)  -  0  (12) 

The  exact  solution  of  this  problem  features  a  shock  layer  that  moves  from  x  ■ 
0  to  x  ■  1  with  unit  speed  and  then  approaches  the  steady  state  solution 
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as  t  ♦  •. 


u(x,t) 


1  -  «-!/* 


The  maximum  error  at  steady  state 

nax|u(j Ax.nAt)  -  Ut' j  |  ,  n  ♦  • 


(13) 


(14) 


computed  by  the  Lax-Weodroff ,  upwind,  LSE,  and  Lax-Friedrlchs  schemes  are 
shown  in  Table  2  for  Ax  •  1/20,  p  ■  6,  a  ■  0.375,  0.75  and  Lx  •  1/20,  p  ■  500, 
a  -  0.475,  0.95*. 

The  centered  difference  scheme  produced  overflows  for  both  p  ■  6  end 
500,  so  no  results  could  be  listed  for  it.  The  Lax-Friedrlchs  scheme  only 
overflowed  for  p  ■  6.  The  LSE  scheme  gives  the  exact  steady  state  solution 
for  this  example  and  the  small  errors  reported  in  Table  2  are  due  to  the 
combined  effects  of  roundoff  and  our  failure  to  reach  the  exact  steady  state. 

These  results  are  very  encouraging;  however,  when  we  examine  the  LSE 
solution  during  the  transient  phase  of  the  solution  the  situation  is  quite 
different  (cf.  Figure  2).  The  LSE  solution  is  overly  diffusive  and  the 
computed  solution  is  not  much  better  than  that  obtained  by  upwind 
differencing.  This  observation  was  also  noted  by  Cresho  and  Lee  [7)  about 
methods  that  are  similar  to  the  LSE  method. 

III.  THE  LINEARIZED  TRANSIENT  EXPONENTIAL  (LTE)  METHOD.  We  would  like 
to  Improve  upon  the  results  of  the  LSE  method  for  transient  problems  and, 
thus,  we  consider  developing  a  method  having  the  form  of  (2)  that  gives  a 
polntwlse  exact  solution  of  (1)  when  f(u)  ■  cu  and  c  is  a  constant.  Since  we 
are  primarily  lntereated  in  obtaining  good  resolution  near  shock  and  boundary 
layers  we  choose  to  solve  (1)  subject  to  Rlemann  initial  data.  To  be 
specific,  for  each  j  and  n  we  compute  as  the  exact  solution  of  the 

initial  value  problem. 

ut  +  cux  -  cux*  ,  t  >  nAt  ,  |x|  <  -  , 

(  UL  »  x  <  (j-l+fi)Ax 

u(x,nAt)  -  (15) 

{  ug  ,  x  >  (j-l+6)Ax 

where  u^  and  ur  are  constants,  A  is  a  constant  on  (0,1)  to  be  determined,  and 
we  assume  that  c  >  0.  We  shall  present  results  for  c  <  0  later. 

The  exact  solution  of  (15)  at  x  -  JAx  and  t  •  (n+l)At  is 

u(jAx,nAt)  -  uR  -  -(ug-UL)erfc  -  *[^"( l-6-o)  (16) 

2  2  T  o 


*A11  numerical  results  were  obtained  in  double  precision  on  an  IBM  4341 
computer  at  the  Benet  Weapons  Laboratory. 
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whereas  the  eolution  of  equetlon  (2)  at  thle  point  la 

..  1 

D^lj  ■  UR  -  -  a(l  -f  tnj_i/2  +  2/p)(ur-ul)  (17) 

The  two  eolution*  will  be  the  eeae  provided  that 

2  1  1  fp 

*n1-l/2  ■  “1  —  +  -  erfc  -%}-  ( 1-6-a)  (18) 

J  pa  2  f  a 

In  this  paper,  we  simplify  (18'/  by  aeswing  that  p/a  »  1  and  approxi¬ 
mating  the  complementary  error  function  by  2H(6+a-l),  where  H  is  a  Heaviside 
function.  Also,  since  snj+i/2  ia  not  determined  by  this  procedure,  we  specify 
it  according  to  equation  (6)  with  coth  p/2  approximated  by  unity.  Thus,  we 
have 


2  1  2 
*nj-l/2  “1  —  +  2(“  H(6+a-l)  -  lj  ,  *nj+l/2  "  1  "  _  <19a> 


When  c  <  0  we  choose 


j-1/2  * 


2 

-1  -  - 

P 


2  1 

znj+l/2  -  -1  -  -  +  2[-  H( 6-a-l)  +  1) 


(19b) 


It  remains  to  specify  6.  One  possibility  is  to  choose  it  to  be  a  random 
variable  uniformly  distributed  on  (0,1),  in  which  case  equations  (2)  and  (19) 
would  yield  a  linearized  random  choice  scheme  [1,2].  A  second  possibility  is 
to  always  select  6-1/2  which  would  give  a  Godunov  [6]  type  scheme.  The 
third  possibility  is  similar  to  a  scheme  suggested  by  Roe  [15]  and  is  the  one 
that  we  have  been  using.  We  begin  by  selecting  6-0;  however,  any  value  of 
6e[0,l)  will  do.  After  each  time  step  we  add  the  magnitude  of  the  Courant 
number  |a|  to  6  and  obtain  a  new  value  of  6.  We  continue  this  process  until  6 
exceeds  unity,  in  which  case  we  replace  6  by  6-1.  The  procedure  has  to  be 
modified  slightly  when  a  is  not  a  constant  and  we  shall  Indicate  how  this  is 
done  shortly;  however,  if  a  is  a  rational  number  of  the  form  p/q  and  c  -  0 
then  equations  (2)  and  (19)  have  the  advantage  of  giving  the  pointwise  exact 
solution  of  the  linearized  Riemann  problem  every  q  time  steps. 

We  refer  to  the  scheme  consisting  of  equations  (2)  and  (19)  and  the  above 
choice  of  6  as  the  linearized  transient  exponential  (LTE)  method  and  we  begin 
by  applying  it  to  the  following  linear  Riemann  problem. 

Example  2: 


ut  +  ux  -  eu*x  ,  t  >  0  ,  |x|  <  - 


u(x,0)  - 


1  ,  x  <  0 
0  ,  x  >  0 


In  this  example,  the  initial  discontinuity  becomes  a  shock  layer  which  travels 
with  unit  speed  in  the  positive  x  direction  while  widening  as  t  increases. 

We  have  computed  the  solution  of  this  problem  by  the  LTE  method  with  Ax  - 
1/20,  p  -  500,  and  a  -  0.75,  0.95.  For  this  value  of  p  and  for  times  less 
than  order  1/c,  the  shock  layer  is  well  contained  within  one  mesh  subinterval 


-  ■  -  '*^t  wjj, 


and  ve  have  plottad  che  locatlona  of  Che  ends  of  this  aublnterval  along  with 
the  exact  poeltlon  of  the  center  of  the  ahock  layer  in  Flgurae  3  and  A  for  o  - 
0.75  and  0.95,  respectively.  Ve  aee  that  the  shock  layer  Is  tracked  exactly 
on  the  average  and  that  ve  obtain  the  polntwlse  exact  solution  every  A  and  20 
tine  steps  when  a  ■  0.75  end  0.95,  respectively. 

For  nonlinear  scalar  probleas  we  stilJ  use  equations  (2)  and  (19); 
however,  ve  now  use  local  values  of  the  Courant  and  Reynolds  numbers  based  on 
a  local  shock  speed.  Thus,  on  each  aublnterval  we  calculate 

anJ-l/2  "  ■nj-i/2At/Ax  t  Pnj-l/2  "  »nj-i/2Ax/e  (21a, b) 

where  fnj-i/2  i»  a  local  shock  speed  which  we  choose  as 

•nj-l/2  "  (fnj  "  fnj-i>/(Dnj  -  0Dj-i)  (22) 

Equations  (21)  are  used  in  equations  (19)  to  calculate  *nj-i/2  and  we  proceed 
as  In  the  linear  case.  After  each  tine  step  we  add 

(mln|anj_i/2 1  +  maxlc^.^l)/! 

to  6  and  obtain  a  new  value  of  6.  Once  6  exceeds  unity  we  again  replace  It  by 

6-1. 


1 


Equation  (22)  gives  the  exact  shock  speed  whenever  U°j  and  Unj-i  satisfy 
the  Ranklne-Bugoniot  Jump  conditions  (cf.  e.g.  Vhlthsm  [18]  and  equation 
(27)).  An  alternate  definition  of  anj-i/2  that  Is  easier  to  use 
computationally,  but  only  gives  the  correct  shock  speed  when  f  Is  at  most  a 
quadratic  function  of  u  is 

*nj-l/2  "  \  +  f'CUVi))  •  (23) 

Example  3:  We  consider  a  Riemann  problem  for  Burgers'  equation 
ut  +  -  (u1 2 * * S)x  -  cu «  ,  t  >  0  ,  |x|  <  •  , 

(  UL  .  *  <  0 

u(x,0)  -  (2A) 

(  uR  ,  x  >  0 

The  exact  solution  of  this  problem  can  be  obtained  by  the  Hopf-Cole  transfor¬ 
mation  and  is  given  In,  e.g.,  Whitham  [18].  Herein,  It  suffices  to  give 
asymptotic  formulas  which  are  valid  for  t/e  »  1.  Thus,  when  u^  >  up  we  have 


1  1  UL  “  UR 

u(x,t)  ~  -(uL  +  ur)  -  -(ul  ”  ur)  tanh  (— - - )(x  -  St)  (25a) 

2  2  Ac 

where 

1 

S  ■  -  (ul  +  ur)  (25b) 


I 


mi 


A 
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end  when  ul  <  ur  we  have 


UL  ,  x/t  <  uL 

u(x,t)  ~  x/t  ,  UL  <  x/t  <  UR  (26) 

UR  »  UR  <  x/t 

Equation  (25)  represents  a  shock  layer  moving  In  the  positive  x  direction  with 
speed  S  and  equation  (26)  represents  an  expansion  fan. 

We  calculated  solutions  with  e  ■  10“**,  Ax  -  1/20,  and  X  ■  0.95  by  the  LSE 
and  LTE  methods  for  a  shock  problem  with  i*l  ■  1,  ur  ■  0  and  an  expansion 
problem  with  ul  ■  -1 ,  ur  -  -1/2.  In  Figure  5  we  compare  the  exact  shock 
position  with  those  calculated  by  the  LSE  and  LTE  methods.  We  define  the 
shock  position  for  the  numerical  methods  as  the  point  where  the  solution  is 
(ul  -  ur)/2  when  linear  interpolation  is  used  to  compute  solution  values 
between  mesh  points.  In  Figures  6a  and  6b  we  plot  the  exact  LSE  and  LTE 
solutions  at  t  ■  0.95  for  the  shock  problem  and  at  t  ■  0.38  for  the  expansion 
problem,  respectively.  The  LTE  method  again  confines  the  shock  layer  to  one 
subinterval  and  gives  the  correct  shock  speed  and  position  on  the  average. 

The  LSE  method  is  overly  diffusive  and  is  giving  the  correct  shock  speed,  but 
the  position  is  wrong  by  about  Ax/2. 

The  situation  is  quite  different  for  the  expansion  fan.  The  LSE  method 
is  still  overly  diffusive;  however,  the  LTE  method  is  representing  the 
expansion  fan  as  a  shock.  This  phenomenon  also  occurs  with  Roe's  scheme  for 
hyperbolic  systems  (cf.  Roe  [16]  and  van  Leer  [17])  and  it  must  be  remedied  if 
these  schemes  are  to  be  useful  on  expansion  problems. 

IV.  SYSTEMS  OF  EQUATIONS.  In  principle  the  LTE  method  consisting  of 
equation  (2),  (19),  and  (21)  may  be  directly  extended  to  vector  systems  of 
the  form  (1)  once  we  have  selected  a  shock  Bpeed  Bnj-i/2»  When  e  -  0  the 
exact  shock  apeed  S  is  determined  by  the  Rankine-Hugonlot  condition 

(ur  -  ul)S  -  f(uR)  -  f(uL)  (27) 

where  ur  and  jjl  are  the  values  of  u(x,t)  on  opposite  sides  of  the  shock.  When 
e  is  nonzero  but  str.all  we  would  like  the  numerical  shock  speed  anj_j/2  "  S 
whenever  the  numerical  solution  iFj.j,  Unj  satisfies  (29). 

In  a  recent  paper  Harten  and  Lax  [8]  suggested  selecting 

•°j-l/2  "  *(fnj  “  fj-lJWj  "  fj-l>  (28a) 

where  t(w)  is  the  linear  functional 

m 

t(w)  -  (V'^j)  -  VWj-jMw  (28b) 

and  V(u)  is  an  entropy  function.  They  show  that  this  choice  gives  unique 
physically  admissible  numerical  solutions  of  their  random  choice  finite 
difference  methods. 
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Roe  (16]  euggeete  ao  alternate  method  of  calculating  >nj-}/2  that  la 
based  on  the  elgenvalut  of  a  matrix  approximating  the  Jacobian  3f/3u. 

*  m 

We  have  not  tried  either  of  these  alternatives,  but  instead  use  the  very 
simple  prescription 

(pj  -  fj-i>T(£°j  -  Tj-i) 

•Vl/2 - - -  (29) 

(pj  -  SQj-i)T(Snj  -  HnJ-i> 

Equation  (29)  gives  the  exact  shock  speed  whenever  and  U°j-i  satisfy  the 
Ranklne-Hugoniot  conditions  (27),  but  it  may  fall  to  give  a  physically 
acceptable  solution. 

Example  4:  We  solve  the  following  impact  problem  for  the  linear  wave 
equation 


ult  "  «2x  ■  0  .  u2t  -  ulx  *  0  , 

uj(x,0)  -  0  ,  U2(x,0)  -  I 


t  >  0  ,  |x|  <  • 

1  ,  x  <  0 


(30) 


{  -i  »  *  >  o 

Here  uj(x,t)  and  U2(x,t)  represent  the  strain  and  velocity  in  two  elastic  rods 
that  impact  each  other  with  unit  speeds  at  x  -  t  -  0. 


We  calculated  the  solution  of  this  problem  by  the  LSE  and  LTE  methods  and 
by  the  EPIC-2  code  [12].  The  latter  is  a  two-dimensional  finite  element  code 
for  elastic-plastic  impact  problems.  Our  results  for  uj  at  t  ■  0.95  obtained 
with  Ax  -  1/10  and  X  -  0.95  are  shown  in  Figure  7.  The  LSE  and  LTE  solutions 
are  typical  of  our  results  on  previous  examples.  The  LTE  method  again  calcu¬ 
lates  the  correct  shock  position  and  speed  with  no  diffusion  or  oscillations. 
The  LSE  solution  is  overly  diffusive,  although  less  so  than  the  EPIC-2  solu¬ 
tion. 


V.  DISCUSSION  OF  RESULTS.  The  LTE  method  appears  to  be  a  very  promising 
acheme  for  shock  problems.  It  le  simpler  to  apply  than  methods  based  on  the 
exact  solution  of  Riemann  problems  [1,2]  and  does  not  suffer  from  the  effects 
of  artificial  diffusion  or  spurious  oscillations.  However,  our  results  are 
very  preliminary  and  there  are  still  many  questions  to  be  answered  and  many 
problems  to  be  overcome.  The  performance  of  the  LTE  method  in  regions  of 
expansion  must  be  improved,  van  Leer  (17]  has  suggested  incorporating 
expansion  fans  in  the  approximate  Riemann  solution  of  Roe's  method  [16],  and 
this  approach  should  work  for  our  LTE  methou  as  well.  Another  possibility  is 
to  base  the  difference  scheme  (2)  on  the  exact  solution  of  (1)  when  f  is  a 
quadratic  function  of  u.  The  solution  of  this  problem  does  contain  expansion 
waves;  however,  extending  this  method  to  systems  of  equations  would  be 
considerably  more  difficult  than  extending  the  LTE  method. 

Both  the  LSE  and  LTE  methods  are  first  order  accurate  lu  time  when  the 
solution  is  smooth,  van  Leer  [17]  has  developed  a  two-step  procedure  that  can 
be  used  to  extend  these  methods  to  second  order  accuracy  and  we  plan  to 
experiment  with  it  shortly. 
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There  le  also  the  possibility  of  developing  implicit  exponentially  fitted 
sod  weighted  schemes,  which  would  be  desirable  when  approaching  a  steady  etate 
and  which  aay  improve  the  phase  characteristics  of  the  LSE  method  (cf.  Cresho 
and  Lee  [7] ) . 

Finally,  we  note  that  the  LSE  and  LTE  methods  can  be  extended  to  higher 
dimensions  by  using  operator  splitting  techniques.  However,  this  aay 
Introduce  some  nuaerical  diffusion. 
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TABLE  1 .  VALDES  OF  i  AND  0  FOR  DIFFERENCE  METHODS 
THAT  HAVE  THE  FORM  OF  EQUATION  (11) 


Method 

a 

0  ■  a(s  +  2/p) 

Centered 

0 

2  a/p 

Lax-Wendrof f 

a 

a2  +  2a/p 

Upwind 

sgn( p) 

a(sgn(p)  +  2/p) 

Linearised  steady 
exponential  (LSE) 

coth(p/2)  -  2/p 

acoth(p/2) 

Lax-Frledrichs 

i/a 

1  +  2a/p 

TABLE  2.  MAXIMUM  ERROR  AT  STEADY  STATE  FOR  EXAMPLE  1. 
AN  *  INDICATES  THAT  THE  COMPUTED  SOLUTION 
PRODUCED  AN  OVERFLOW. 


Method 

p  -  6 

p  “  500 

a  -  0.375 

0.75 

a  -  0.475 

0.95 

Lax-Wendrof £ 

2.9  E-l 

1.6  E-3 

3.5  E-l 

2.3  E-2 

Upwind 

2.0  E-2 

2.0  E-2 

2.0  E-3 

2.0  E-3 

Linearised  steady 

2.5  E-l 4 

2.5  E-l 4 

8.5  E-8 

2.1  E-12 

exponential  (LSE) 

Lax-Frledrichs 

* 

* 

3.6  E-l 

2.8  E-2 
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Figure  1.  legion  of  linear  stability  for  equation  (11)  and  curve*  of  0  vs.  a 
for  the  centered  difference  (C),  Lax-Wendroff  (LW),  upwind 
difference  (U),  linearised  steady  exponential  (LSE),  and 
Lax-Friedrlchs  (LF)  methods. 


Figure  2.  Ooaparlson  of  exact  and  LSE  solutions  of  Example  1  at  t  •  0.475. 

Calculations  were  performed  with  Ax  -  1/20,  a  -  0.95,  and  p  -  500. 


I  ’ 
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Figure  3.  Exact  shock  layer  position  and  the  location  of  the  subinterval 

containing  the  shock  layer  calculated  by  the  LTE  method  for  Example 
2  with  Ax  ■  1/20,  a  -  0.75,  and  p  •  500. 
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Figure  5.  Exact  shock  leyer  position  end  those  calculated  by  the  LSE  and  LTE 
methods  for  Example  3  with  c  -  10"\  Ax  -  1/20,  and  a  -  0.95. 
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Figure  6.  Exact,  LSE,  and  LTE  solutions  of  Example  3  with  c  •  10*“  ,  Ax  • 

1/20,  and  o  -  0.9S.  In  (a)  we  show  the  solution  it  t  ■  0.95  of  a 
ahock  problem  with  uj,  ■  1,  ur  -  0,  and  in  (b)  we  show  the  solution 
at  t  ■  0.38  of  an  expansion  problem  with  ur  ■  -1,  ur  •  -1/2. 
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Figure  7 
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Oonpariaoa  of  exact,  LSE,  LTE,  and  EPIC- 2  aolutlona  for  uj  of 
Exaaple  4  at  t  •  0.95  with  Ax  -  0.1  and  At  •  0.095. 
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